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Abstract 

Continuum solvent models have become a standard technique in the context of electronic struc- 
ture calculations, yet, no implementations have been reported capable to perform molecular dynam- 
ics at solid-liquid interfaces. We propose here such a continuum approach in a DFT framework, 
using plane-waves basis sets and periodic boundary conditions. Our work stems from a recent 
model designed for Car-Parrinello simulations of quantum solutes in a dielectric medium [J. Chem. 
Phys. 124, 74103 (2006)], for which the permittivity of the solvent is defined as a function of the 
electronic density of the solute. This strategy turns out to be inadequate for systems extended in 
two dimensions: the dependence of the dielectric function on the electronic density introduces a 
new term in the Kohn-Sham potential which becomes unphysically large at the interfacial region, 
seriously affecting the convergence of the self-consistent calculations. If the dielectric medium is 
properly redefined as a function of the atomic coordinates, a good convergence is obtained and the 
constant of motion is conserved during the molecular dynamics simulations. Moreover, a significant 
gain in efficiency can be achieved if the simulation box is partitioned in two, solving the Poisson 
problem separately for the "dry" region using fast Fourier transforms, and for the solvated or "wet" 
region using a multigrid method. Eventually both solutions are combined in a self-consistent pro- 
cedure, and in this way Car-Parrinello molecular dynamics simulations of solid-liquid interfaces can 
be performed at a very moderate computational cost. This scheme is employed to investigate the 
acid-base equilibrium at the Ti02-water interface. The aqueous behavior of titania surfaces has 
stimulated a large amount of experimental research, but many open questions remain concerning 
the molecular mechanisms determining the chemistry of the interface. Here we make an attempt 
to answer some of them, putting to the test our continuum model. 
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I. INTRODUCTION 



The structure and the reactivity of sohd surfaces have become major subjects of study 
in chemistry and condensed matter physics, and are at the core of much of the research 
conducted in materials science. In this context, electronic structure methods, in particular 
density functional theory (DFT) calculations, have played a fundamental role in the inter- 
pretation and the validation of the data coming from the different surface spectroscopies 
and microscopies, often providing new insights on top of those detaching from the available 
experimental techniques.^ Throughout the past two decades a large number of DFT sim- 
ulations has been reported on metallic and semiconducting surfaces, contributing precious 
information concerning atomic and electronic structure, thermodynamics, and reactivity.-"- 
With very few exceptions, these simulations considered a slab in the gas phase. For a broad 
range of applications, however, the relevant phenomena occur in the presence of a liquid 
phase, as is often the case in processes related to electrochemistry and catalysis. The real- 
ization of liquid phase DFT simulations is therefore a much pursued objective — especially 
when many of the standard X-ray techniques like XPS or SAXS are unsuited to provide 
atomic scale information in solution — , but the inclusion of the solvent considerably in- 
creases the cost of first-principles calculations of periodic surfaces, and is therefore a rather 
uncommon practice.- 

In an explicit solvation approach, the vacuum space in the simulation box is filled with 
solvent molecules. The subsequent growth in the size of the system is in part responsible for 
the increased computational expense, but, more importantly, there is the fact that a static 
picture is a very poor representation of the hquid state. Any solute admits a huge number 
of possible configurations for the solvent molecules around it, associated with multiple local 
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minima, the net solvation effect arising from a weighted average of all these.- A geometry 
optimization would lead to a solid or glassy phase corresponding to one of these minima in 
the multidimensional potential energy surface, resulting in a dielectric screening typically 
different from the static limit observed in the liquid state. Hence, to capture the solva- 
tion effect, it is necessary to perform extensive statistical sampling involving either lengthy 
molecular dynamics or Monte-Carlo simulations. Alternatively, it is possible to resort to 
the so called continuum (or implicit solvent) models, in which the solvent molecules are 
replaced by a dielectric medium surrounding the solute and exhibiting the static screening 
of the solution.-i^ In this way, the polarization induced by the solvent is introduced in an 
averaged fashion, and the cost of the computation gets closer to the corresponding cost in 
vacuum by drastically reducing the number of degrees of freedom. On the other hand, the 
representation of the solvent structure is omitted, disregarding any possible solute-solvent 
specific interactions. Still, to overcome this problem, all or part of the first solvation shells 
can be included explicitly, with the dielectric medium extending beyond the limits of this 
cluster comprising the solute plus a few solvent molecules. In any case, the continuum model 
has a long tradition in quantum chemistry and has proved reliable and efficient to extract 
properties in solution of a large variety of molecular systems.-"- 

In recent years a small number of implementations of the continuum model has been 
proposed in the context of density functional theory, plane-waves basis sets, and the Car- 
Parrinello method.—"— To the best of our knowledge, none of them have been employed 
in molecular dynamics simulations of periodic surfaces.— Here we revise the electrostatic 
continuum model of Fattebert and Gygi,— proposing a new definition for the dielectric 
function to make the model suitable for the treatment of periodic slabs. The original for- 
mulation, in which the permittivity of the solvent is determined by the charge density of 
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the solute, has proved successful in the simulation of molecular and ionic solutes,— and 
of extended systems like polymers, with periodicity in one dimension.J^ In the case of solid 
surfaces, however, the convergence to the electronic ground state turns out to be impaired. 
In the sections that follow, we attribute this failure to a term in the Kohn-Sham potential 
originating in the dependence of the dielectric constant on the electronic density, and, to 
circumvent this issue, we redefine the permittivity as a function of the atomic positions. In 
this way the smooth electronic convergence is restored, along with a good conservation of 
the total energy in the molecular dynamics runs. The electrostatic problem in the presence 
of the dielectric medium is efficiently addressed with a multigrid method,—"— which solves 
the Poisson equation in real space. Moreover, a significant increase in efficiency can be 
obtained if the simulation box is partitioned in two, solving the Poisson problem separately 
for the "dry" region using fast Fourier transforms, while restricting the multigrid treatment 
to the solvated or "wet" part of the supercell (and then combining both solutions in a self- 
consistent fashion). Using this scheme it is possible to carry out Car-Parrinello molecular 
dynamics simulations of solid-liquid interfaces at a computational cost exceeding by just 
a small factor the one corresponding to vacuum. We employ this approach to investigate 
the proton exchange at the anatase-water interface, which is a key process in the acid-base 
equilibrium of Ti02 surfaces. The ubiquity of titania in solid-liquid applications has incited 
the emergence of empirical models to assess the protonation and the dissociation of terminal 
groups on the Ti02 surface in aqueous environments.—!^ These simple models have been 
broadly used among experimentalists for the interpretation of their data, but a molecular 
level description that accounts for the effect of structure has not been established. In this 
context, the present simulations intend to provide a first-principles glance of the microscopic 
mechanisms governing the chemistry of this interface. 
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II. MODEL AND METHODOLOGICAL DISCUSSION 



A. About the present implementation and the computational parameters 

This model has been implemented in the public domain Car-Parrinello parallel code in- 
cluded in the Quantum-ESPRESSO package,— based on density-functional theory (DFT), 
periodic-boundary conditions, plane-waves basis sets, and pseudopotentials to represent the 
ion-electron interactions. All calculations reported in this work, unless otherwise noted, 
have been performed using the PW91 exchange-correlation functional^^ in combination with 
Vanderbilt ultrasoft pseudopotentials.-^^ The Kohn-Sham orbitals and charge density were 
expanded in plane waves up to a kinetic energy cutoff of 25 and 200 Ry respectively. Peri- 
odic slabs of four layers width representing the (101) surface of the anatase structure were 
computed using gamma-point sampling in supercells of size 10.24 A x 7.56 A x 17.82 A. 
For finite temperature (not damped) Car-Parrinello molecular dynamics simulations, an 
electronic mass of 400 a.u. and a time step of 0.17 fs were adopted. During geometry relax- 
ations and dynamics, all atoms were allowed to move, with the exception of those belonging 
to the two inner layers, which were fixed in their bulk positions. 

B. The continuum solvation model 

Within the continuum approach, the solvent is represented as a dielectric medium sur- 
rounding a quantum-mechanical solute confined in a cavity. In particular, we consider the 
polarizable continuum model, in which the dielectric medium and the electronic density re- 
spond to the field of each other in a self-consistent fashion. This interaction provides the 
electrostatic part of the solvation free energy, AG^i, which is the dominant contribution 
for polar and charged solutes. The cavitation energy AGcav is defined as the work in- 
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volved in creating the appropriate cavity inside the solution in the absence of solute-solvent 
interactions. Electrostatic, cavitation, and dispersion-repulsion effects are modeled as sep- 
arate contributions,- and the free energy of solvation is regarded as the sum of these three 
terms {AGsoi = AG^i + AGcati + ^Gdis-rep)- Thermal and pressure dependent terms can 
also be included but are usually negligible. We note that this decomposition is inherent to 
the model, being AGsoi the only measurable quantity. 

C. Previous implementation 

The starting point for this work is the implementation reported in reference ~: in the 
following paragraphs we revisit those aspects of the preceding version that are essential to 
the present development. First, we note AGei and AGcav are considered explicitly, while 
AGdis-rep, less relevant for solutes of moderate size, is largely seized as part of the elec- 
trostatic term by virtue of the parametrization. The electrostatic interaction between the 
dielectric and the solute is calculated as proposed by Fattebert and Gygi,^^-'-^^ who define the 
permittivity e of the solvent as a function of the electronic density p. Within a pseudopo- 
tential framework and in periodic boundary conditions, the Kohn-Sham energy functional^^ 
can be written as 

E[p] = T[p] + E,,[p] + <f){r)ptot{r)dr 

where T[p] corresponds to the kinetic energy of the electrons and E^dp] to the exchange- 
correlation energy. The last four terms on the right hand side account for the total electro- 
static energy in a periodic crystal, gathering all Coulombic interactions involving electrons 
and nuclei, only omitting for simplicity the non-local part of the pseudopotential (for a 
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detailed derivation see reference^ or the Appendix of reference^). The electrostatic formu- 
lation in Eq. (1) arises from the Ewald sum of point charges, which requires to introduce the 
ionic densities p/(r — Ri) consisting of Gaussian distributions of negative sign that integrate 



to Z/, the total charge of the pseudo-ion: p/(r — Ri) = exp^— ^''(j^^l ), with Rj 

the width of the Gaussian associated with the site /. 

The third term on the right is conventionally called the Hartree energy: 



Eh = ^J (t>{r)ptot{r)dr 



where ptot is just the sum of the electronic plus the ionic densities, ptotij) = p(i") + 
Y^i Piij^ ~ R'l)) while the electrostatic potential 0[p] is the solution to the Poisson equation 

VV(r) = -A^ptot{v) . (2) 

In the presence of a dielectric continuum with a permittivity e[p], the Poisson equation 
becomes 

V-(e[p]V0(r)) = -47rp,,i(r) . (3) 
Using Eq. (3), the formula for the Hartree energy Eh-, can be integrated by parts to yield: 

EH = l:^je[p]{Vct>{v)fdv. (4) 

The functional derivative of Eh with respect to the charge density is added to the other 
contributions of the energy (the exchange-correlation, the local and non-local parts of the 
pseudopotentials) to set up the Kohn-Sham potential V^^'^fp]. 

8Eh, 



6p 



r 



(r) + K(r), (5) 



K(r) = (V0(r))2|i(r). (6) 

OTT op 



We provide the derivation of Eq. (5) and (6) in the Appendix A, since it is not given in 
any of the previous references. The dielectric medium and the electronic density respond 
self-consistently to each other through the dependence of e on p and vice versa. 

In quantum chemistry continuum models as PCM,-i2I the dielectric constant e is taken to 
be 1 inside the cavity, and a fixed value outside. For molecular dynamics applications, such 
a discontinuity needs to be removed to calculate accurately the analytic derivatives of the 
potential with respect to the ionic positions. Besides, in the particular case of plane-waves 
implementations based on real space grids, a smoothly varying dielectric function is more 
appropriate for numerical reasons. Fattebert and Gygi proposed the following smoothed 
step function for the dielectric: 



This function asymptotically approaches (the permittivity of the bulk solvent) in regions 
of space where the electron density is low, and 1 in those regions where it is high (outside 
the solvation cavity). The parameter po is the density threshold determining the cavity size, 
whereas (3 modulates the smoothness of the transition from eoo to 1. 

On the other hand, the cavitation energy is computed separately as the product between 
the surface tension of the solvent and the area of the cavity. In this implementation this 
term remains unchanged; details about the calculation of AGcav can be found in reference — . 

D. Calculation of periodic slabs: the problem in the convergence 

When the scheme described above is applied on structures extended in two dimensions, 
it typically fails to achieve self-consistency. We found the electrons heat up during the 
Car-Parrinello relaxation of the wavefunction, causing the total energy to diverge, as shown 
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in Fig. [H In other cases the kinetic energy of the electrons is observed to decrease at 
the beginning, but thereafter experiences irregular oscillations without ever reaching zero. 
Even if convergence is enforced through the use of a stringent algorithmic strategy, the 
solution obtained is not reliable, e.g., it depends on the starting conditions or the damping 
parameters. 

The source of this erratic behavior can be tracked to the inclusion of the term K in the 
Kohn-Sham potential. This term, defined in Eq. (6), arises from the dependence of the 
permittivity on the electron density. Fig. [T] illustrates the impact of this contribution on 
the convergence of a Car-Parrinello electronic relaxation, by showing that an acceptable 
dynamics is recovered if is neglected.— In particular, these results correspond to the 
anatase-water structure described in section III, but the convergence of other systems is 
degraded in a similar way. The reason for the instability associated with can be somehow 
appreciated in the two dimensional plot of Fig. [21 which displays the variation of in the 
2;-direction (the one perpendicular to the surface) at fixed x and y. These coordinates have 
been chosen in coincide with the position of a Ti atom of the surface. The same graph shows 
the total Kohn-Sham potential, excluding V^, so that the relative contribution of this term 
can be clearly seen. The values are exceptionally large at the solid-liquid boundary, even 
larger than any of the other contributions to the effective potential. Referenced on the right 
axis, the dielectric function e[p] is also shown. Since depends on de/dp (Eq. 6), the peaks 
occur at the region where the charge density decays abruptly, producing a rapid variation 
in e[p] from the value in the solid to the value in the solvent. Adopting the model function 
of Eq. (7) we have 

dp^' Po (l + (p(r)/poW 
This expression goes to zero when p(r) ^ po or p(r) ^ po, and is dominated by l/po 
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otherwise.— Then, the extreme values of K at the sohd-hquid interface ultimately originate 
in the magnitude of the density threshold po- Unfortunately, this parameter may not be 
freely tuned, but is set in combination with (3 to fit the experimental data, and no reasonable 
choice of (po, /?) can be made as to prevent the blowup of V^. Alternative model functions to 
represent the dielectric were considered, e.g., with e[p] exhibiting a Gaussian or trigonometric 
decay with p. None of these functions entailed any significant improvement, as far as all of 
them have in common a sudden change associated with the transition from 1 to eoo, which 
redounds in large values of de/dp at the interface. The transition can be smoothed enough 
as to avoid the sharp peaks in V^, but in doing so the solvation effect is ruined. Clearly, the 
discussed behavior is a consequence of the dependence of e on p, irrespective of the kind of 
function chosen to model the dielectric. It should be noticed that the inclusion of K does not 
appear to disturb the convergence or the energy conservation in the case of finite systems, 
neither in the case of polymers (extended in one dimension).— It is seemingly because of the 
bi-dimensional symmetry of periodic surfaces that K turns out to spoil the Car-Parrinello 
dynamics, by propagating the observed perturbation of the Kohn-Sham potential throughout 
a full, extended plane. 

E. A position dependent formulation 

It is fair to wonder about the physical meaning of this strong perturbation of the potential 
contained in V^. In the absence of a dielectric, or for a dielectric defined independently of 
the charge density, the functional derivative of Eh with respect to p turns out to be equal to 
the electrostatic potential 0[p] (see reference^). The additional term emerges from the 
dependence of Eu on e, whenever e is modeled as a function of the charge density. Since 
de/dp < (in Eq. (8), eoo > 1), K is of a repulsive character. The choice of p to dehmitate 
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the region filled by the continuum induces a response of the potential to counteract the 
strong discontinuity imposed by the same p. is such a response, which, throughout the 
self-consistent procedure, opposes the perturbative drift enforced by the model. 

Thus, the instability of the model emanates from the use of the self- consistent charge 
density to define the dielectric. A possible alternative would be to use instead a non-self 
consistent, or "fake" density, since the role of the charge density in this context is simply to 
shape the dielectric medium. As will be discussed below, such an option can be equivalent 
to define e as a function of the atomic coordinates, which is the usual strategy in quantum 
chemistry methods. The idea of a dielectric determined by the atomic positions may be less 
attractive from a physical viewpoint: the size and shape of the cavity depend on the identity 
of the atoms only, and is not modulated by the electronic structure or the environment; the 
polarization of the solvent is lost, and, on top of these, many more parameters are needed — 
at least one per atom. In practice, however, it is possible to choose a dielectric function based 
on the molecular coordinates that closely reproduces the p-dependent solvation, because the 
effect of the self-consistency and the polarization of the solvent on AG^o/ is quite minor. 

In this way, we keep the expression for e given in Eq. (7), but feed it with a dummy 
density 7(r) determined by the ionic positions Rj, 



where Rl^w van der Waals radius for the corresponding species. Hence the dielectric 

function takes the following form: 



Using this definition, the transition of e(7(r)) between 1 and etx, is centered around the 
van der Waals radius. Aside from -R^^^, which values are tabulated, (3 is left as the only 
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adjustable parameter to fit the experimental solvation energies. Fig. [3] shows the aspect of 
the dielectric function around an oxygen atom for different (3. This parameter must be large 
enough to ensure most of the transition occurs within a small length window, but at the 
same time the numerical accuracy needs to be preserved, so there is an upper bound for j3 
which depends on the given grid size. 

Within this framework, the electrostatic contribution to l^^'^[p] is simply the electrostatic 
potential and the stability of Car-Parrinello dynamics in periodic slabs is recovered. 
Table I presents, for several neutral and charged solutes, a comparison between the values 
of AG^o; obtained with the dielectric functions of Eq. (7) and Eq. (10), respectively. The 
small differences proceed exclusively from AGe/, since ^Gcav is the same in both cases 
{/S.Gsoi = ^Gei + AGcat))- As mentioned above, with a proper choice of [3 the position 
dependent dielectric is able to provide solvation energies in close agreement with the previous 
model and with experiments. The results shown correspond to (3 = 2.4. 

The explicit dependence of e on the ionic positions involves a new contribution to the 
forces arising from the derivative of Eh with respect to Ri, which must be taken into account 
to perform conservative molecular dynamics simulations. After some manipulation involving 
the use of Eqs. (3) and (4), this derivative can be expressed as follows (the full derivation 
is given in Appendix B): 



The computation of the first term on the right is straightforward since we know, from Eqs. 
(9) and (10), the explicit dependence of e on Ri: 

with r a generic coordinate x, or Ri = (xq, 2/o, -^o) and R= \y — Ri|. 




(11) 
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On the other hand, if (j){r) is transformed to Fourier space so that (j){r) = J2g 0(G)e**^"" 
(see next section), the second term in Eq. (11) can be evaluated as 

/ '^^")^9ir'^' ^ ^ (G)p,(G)e-«^' (13) 

with p/(G) the form factor of the ionic densities, p/(r — Rj) = J2g P/(G)e~**^'"e~**^^'. 

To ensure the conservation of the total energy during the molecular dynamics simulations, 
the contributions given in Eqs. (12) and (13) must take the place of the derivative of the 
Hartree energy in the absence of the dielectric: 

||._4.nE>G(^)MG)e— . (14) 

F. The multigrid scheme and a mixed strategy to solve the Poisson problem 

In standard plane waves codes based on real space grids, the electrostatic potential 0(r) is 
obtained from the Poisson equation (2), which can be efficiently inverted with the use of fast 
Fourier transforms (FFT). Both the total charge density ptoti^) and 0(r) can be expanded 
in plane waves, 

Ptotir) = Y: p(G)e^^^ 0(r) = 0(G)e^«^ 

G G 

Replacing into the Poisson equation = —Anptot, and equating coefficients: 

Att 

GV(G) = 4vrp(G) , 0(r) = ^ ^T^pW"^'' (15) 

G ^ 

Unfortunately, the Poisson equation in the presence of an arbitrary dielectric, Eq. (3), can 
not be addressed in the same way, and an alternative numerical scheme is required. To that 
end, we have implemented from scratch a sixth-order multigrid code,^'^ which solves in real 
space the Poisson equation with non-constant coefficients and periodic boundary conditions. 
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Eq. (3) can be rewritten as: 

dt d(j) ^ dt d(j) ^ dt d(j) ^ / d'^e ^ d'^e ^ 9^e\ ^ 
dxdx dy dy dz dz \dx'^ dy'^ dz'^ J 

This equation is developed in finite differences, expanding the derivatives of and e to 
sixth-order according to the following relations for the gradient and the Laplacian: 

^ = ^ E «n ^i+n + 0{h') (17) 

8'^ f(r) 1 ^ 

^^-J^2i:/nU.^r. + 0{h') (18) 

where r is a generic coordinate x, y, or z; h is the grid spacing in the direction r; and u is 
the discretization of a continuous function /(r), representing 0(r), p(r), or e(r). Ui refers 
to u evaluated at a mesh point associated with r, while corresponds to a neighboring 
point n positions to the right in the direction r (if n < 0,Uj+„ is located to the left of Ui). 
The coefficients and are given by: 

3 3 1 

ao = 0, ai = -, a2 = - — , as = — , "-n = -^n 

4 zU bU 

49 27 27 2 

^°^"18' '^^=18' ^^^"180' ^^=180' = 

In the case of e(r) = 1 for all r, this method provides a solution for 0(r) which is indistin- 
guishable from the one obtained with FFT. It also demonstrated an excellent performance 
when tried on functions with non-constant coefficients and oscillation frequencies compa- 
rable to those of interest. For example, for = e""*" and e = e"'"' (0.5 < a,b < 2.0) the 
relative error in 0(r) was less than 10~^ in a mesh of 80x80x80 points. 

At the initial steps of a molecular dynamics simulation, the convergence of the potential 
may require 15-30 multigrid cycles. Given the self-consistent nature of the procedure, how- 
ever, after a few time steps the number of cycles necessary to reach convergence is typically 

15 



decreased to less than five. Even so, the multigrid algorithm is still significantly more expen- 
sive than FFTs. Propitiously, multigrid methods can be adapted to any kind of boundary 
conditions, and this feature can be exploited to reduce the size of the mesh involved. To 
implement this idea, a region in the supercell — preferably the slab — must remain inacces- 
sible to the solvent, so that e(r) = 1 within it. In practice, the dielectric function of Eqs. 
(7) or (10) is not diffuse enough as to encompass all the volume of the slab — if it were, the 
solvation effect would fade at the molecular boundaries — , so the solvent occupies the inter- 
stitial space left by the atomic structure (Fig. H^). This is inconvenient not only because 
it increases the numerical complexity of the problem, but also because it is not physical, 
i.e., the solvent does not penetrate the atomic structure of the surface. A simple device to 
exclude the solvent from the solid interspaces is to modify 7(r) in the following way: 



where zj is the z-component of Rj (we recall z is the coordinate perpendicular to the surface; 
it is zero at the bottom of the unit cell and maximum at the top). The factor ^Ef^ produces 
a rapid increase in 7(r) underneath the atom /, which saturates the value of e(r) for all atoms 
/ located below zum- Thus, only the upper face of the slab is in contact with the solution, 
the dielectric function becoming equal to 1 throughout the lower section of the supercell. 
Fig. IDd depicts e(r) when is chosen equal to the z coordinate of an ion belonging to the 
second layer. 

Under these conditions, it is possible to solve the Poisson problem in two steps: first, 
the potential 0(r) is computed for e = 1 in the full domain VL using FFTs; secondly, 0(r) 
is recalculated in the presence of the dielectric with the multigrid method but in a smaller 
domain fii, imposing the boundary conditions proceeding from the FFT solution. The sub- 




(19) 
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domain Qi must be chosen in such a way that it overlaps with a sub-space where e(r) = 1, 
so that the boundary conditions for the multigrid are given by the solution of the Poisson 
equation in vacuum. Then, the total potential is constructed from the combination of the 
multigrid solution in fli plus the FFT solution in the rest of the space. The full procedure 
is summarized in Scheme 1. 



Q 



FFT 
£=1 



0(0) 




MG ^ 

— rT0 (Q) 

£=£(r) 







MG 







FFT 



Scheme 1 

In proposing this approach, we assume the dielectric effect is mostly local and does not 
perturb the electrostatic potential deep inside the slab. More formally, this implies that 
even in the presence of the dielectric medium, there is a region in the real space grid where 
the self-consistent potential is exactly the same as the potential corresponding to vacuum. 
It is therefore important that Jli extends sufficiently down into the surface as to reach that 
region. With this mixed scheme, the size of the mesh to be submitted to the multigrid 
routine can be reduced up to 40 %. 

III. APPLICATIONS: MOLECULAR DYNAMICS AT THE TIO2 INTERFACE 

It is well known that surface groups of most inorganic oxides ionize in solution, exhibiting 
the following equilibria: 
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+ H3O+ ^ M-OH^+^ + H2O ^ + HO" 

with z the surface group charge, which can be negative, positive, or zero, depending on the 
nature of the oxide.~ Understanding the acid-base behavior resulting from these equihbria 
is crucial in almost every application of these materials in solution. Isoelectric points of 
many oxides have been known for years, however, it is very difficult to probe the surface 
of bulk materials or nanoparticles in solution, and most of the data collected corresponds to 
the average behavior of the different surfaces exposed in a given experiment. More recently, 
researchers have sought to take advantage of density functional theory to establish the 
degree of dissociation and protonation at different titania-water interfaces with an explicit 
representation of the solvent.—"— For the reasons already discussed, such an approach is 
costly and has been employed only in a limited number of cases. In what follows, we apply 
our continuum solvent model to characterize the hydrated (101) surface of the anatase 
structure of Ti02, which is possibly the most stable and abundant.- The adsorption of 
H2O on perfect Ti02 surfaces in the gas phase has been thoroughly investigated using 
both experimental and theoretical approaches.—"— In the case of the (101) face of anatase, 
there is consensus in the fact that molecular adsorption of water is thermodynamically the 
most stable. Electronic structure calculations suggest that the difference between the two 
possible adsorption modes — molecular versus dissociated — is of nearly 10 kcal/mol.™ We 
have performed calculations in four layers slabs representing the (101) surface of the anatase 
structure. As previously reported, our own calculations in vacuum summarized in Table 
II show that at different water coverages, the molecular pathway is the most favored. The 
same trend prevails in the presence of the solvent, although the interaction energies with 
the surface turn out to be significantly lower. Table II presents the results from geometry 
optimizations of the fully hydrated surface embedded in the continuum dielectric. The 
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observed weakening of the interaction with the oxide with respect to vacuum is a consequence 
of the stabihzation of the H2O molecules in the polar environment, and can be understood 
in terms of a competition between the bulk solvent and the surface for the water molecules. 
Despite this lower affinity, the massive presence of water from the liquid phase will displace 
the equilibrium toward the formation of an adsorbed monolayer. The energy difference 
between the two kinds of mechanisms remains about the same as in the gas phase, the 
dissociative adsorption becoming exothermic. In solution, however, dissociation is likely to 
occur, controlled by the pH of the medium (see below). 

The quantitative effect of pH on the ionization of the surface is quite difficult to assess 
from first principles simulations, since a huge supercell would be needed to have a meaningful 
representation of the proton concentration in the system. In this preliminary study, we limit 
ourselves to examine the proton exchange between an adsorbed water molecule and an 
hydroxyl anion from the solution, using molecular dynamics at 300 K. This computational 
experiment is meant to provide a qualitative picture of the abstraction of a proton from 
the surface in the presence of OH", illustrating at the same time the performance of the 
continuum solvent method. The inset of Fig. O shows the evolution of the dynamics through 
a sequence of photos, starting at an initial configuration in which an 0H~ group exhibits 
an H-bond with an adsorbed water molecule. Early in the simulation, the covalent 0-H 
bond in H2O is disrupted, to leave an hydroxyl function on the surface and a newly formed 
water molecule that soon wanders around the liquid phase. Fig. \5\ presents the Oq-Hq and 
Ofe-Ha distances, where 0^ and 0^ are the oxygen atoms of the (initially) absorbed water 
molecule and the hydroxide, respectively, and is the abstracted proton. For this reaction 
to occur in a reasonable time as to be within reach of molecular dynamics simulations, the 
starting geometry must be chosen appropriately. As a matter of fact, had the simulation 
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been started from a random configuration, the 0H~ anion might have explored the supercell 
for several picoseconds without ever reacting with the water molecule. In the absence of the 
solvent, instead, the unscreened interaction between the hydroxide and the surface leads to 
an immediate reaction. This distinctive behavior is displayed in Fig. [HI which presents the 
Oft-Hfl distances for two simulations, one in vacuum and the other in solution, started from 
the same geometry and with identical computational parameters. The observed contrast 
between the two dynamics evinces how the dielectric medium stabilizes the hydroxide in the 
liquid phase. The use of this kind of methodology in combination with weighted importance 
sampling techniques (e.g.. Umbrella sampling)^ could provide estimates for the acid-base 
equilibrium constants corresponding to different oxide surfaces and phases. We believe this 
direction, even though beyond the scope of the present work — with an emphasis on the 
methodological conception — aims to a very appealing ground for future research. 



IV. CLOSING REMARKS 

We have shown that a dielectric medium defined as a function of the self-consistent charge 
density provokes a strong response in the effective potential, which in solid-liquid systems 
may spoil the convergence of the Car-Parrinello electronic dynamics. Such a response can 
be avoided with a dielectric based on a non self-consistent charge, preserving in this way the 
potential and allowing for conservative molecular dynamics simulations. This approach is 
equivalent to have a position-dependent permittivity, and therefore a new term in the ionic 
forces must be considered. 

The methodology presented here is a powerful instrument to explore the thermodynamics 
and the reactivity of surfaces and nanoparticles in solution. Replacement of the solvent 
molecules by a dielectric continuum may neglect the structural features of the liquid phase, 
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but it does capture the essential polarization effect of the medium. This is manifest, for 
instance, in the charge of the hydroxyl group: if an additional electron is added to a neutral 
system consisting of a slab plus a (distant) OH moiety in vacuum, a significant portion of 
the charge flows to the solid phase. Noticeably, in the presence of the solvent, the excess 
electron spontaneously localizes on the hydroxyl. This is quantitatively depicted in Fig. [3 
A compelling application for this continuum solvent scheme, as well as a natural contin- 
uation of this work, would be the characterization of the adsorption energies and geometries 
of water and other species on the different surfaces of titanium dioxide. We deem especially 
worthwhile the calculation of the reaction energies for the kind of equilibria mentioned above, 
e.g., Ti-0H2 + 0H~ Ti-OH~ + H2O, as a function of the surface structure. Work in this 
direction is now in progress. At this point it should be noted that it would be not feasible 
to have an estimate of these quantities without the solvation model: in the absence of the 
dielectric, the interaction between the surface and the 0H~ ion (or between the charged 
slab and the water molecule) is extremely dependent on the distance separating them, and 
therefore it is not possible to establish unequivocally the energies for reactants and prod- 
ucts. When the dielectric is included, the long range interaction between charged and polar 
fragments is efficiently screened, and the total energy of the system becomes independent 
of the position of the molecule (or the ion) with respect to the slab. This property makes 
possible to evaluate reaction energies on periodic surfaces in solution which could not be 
calculated by other means, except perhaps with extensive molecular dynamics simulations. 
Aside from these, the scheme presented here would be particularly useful to assess the role 
of the solvent in a great diversity of problems in surface chemistry, including the effects on 
structure, on vibrational frequencies, or on charge transfer phenomena, among many others. 
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VI. APPENDIX A: FUNCTIONAL DERIVATIVE OF V^ 

Following Parr and Yang,~ the functional derivative of Eh[p\ with respect to p(r), which 
we will denote is defined by the relation 



lim 



EH{p + Xf)-EH{p) 



j'-l^n^y. (20) 



A 

where /(r) is arbitrary. From the definition of Eh[p] we can write 

2Eh{p + A/) - 2Eh{p) = 1 0p+A/(p + A/) dv - J 0,prfr 

= J (pp+xfpdv + X j (pp+xffdr - J 0pprfr = A J (l)p+xf f dr + J {(pp+xf - <Pp) P dr. (21) 

In the last term above p can be rewritten using Eq. (3), and the resulting expression can 
be integrated by parts (from now on we omit dr from the integrand for conciseness): 

j {^p+xf - (Pp) P = -j:^ J {(Pp+xf - <Pp) V (epV0p) 

^ 'L j ^'t>p+xf^p^<Pp - ^/ V0pepV0p 
4^ / V0P+A/ (ep - ep+xf) V0p + ^ / ^<Pp+xf ^p+xf^^p - ^ / ^^p ^p^^Pp 



^ / V0P+A/ (ep - ep+A/) V0P + /(P + A/)0p - ^ / ep(V0, 



2 

PJ ■ 
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To arrive to the last expression, ep+A/ was added and subtracted in the third equahty, and 
then the second term was integrated by parts. This resuh can now be inserted in Eq. (21): 

2Eh{p + A/) - 2Eh{p) = xJ (l^p+xf f + ^J V0p+A/ (cp - Cp+A/) V0p + A | /^^ 
Dividing Eq. (22) by 2A, and taking the hmit (A — > 0), we get to the final outcome: 



lim 



EHip + Xf)-EH{p) 
A 







= lim 




A^O 





(cp+A/ - 



By comparison with Eq. (20), it is easy to see that — 4>p — ^C^^p)'^^- 

VII. APPENDIX B: DERIVATIVE OF Eh WITH RESPECT TO THE IONIC PO- 
SITIONS 

In the absence of a dielectric (e = 1), Eq. (15) can be inserted in the expression for the 
Hartree energy to give 

E„ = 2na^^^ (24) 

where p(G) are the Fourier coefficients for the expansion of ptot{r), p(G) = Pe(G) + 
J2i p/(G)e~*^^^ Since the only explicit dependence of p(G) on {Ri} is through the struc- 
ture factor (e^*^^'), the derivative of Eq. (24) with respect to the atomic positions is 
just 

^ = _4..S'g(^)mG).— . (25, 

In the presence of a dielectric medium determined by the ionic coordinates, Eq. (24) does 
not hold. To calculate dEn/dKi we replace e[p] for e(Ri) in Eq. (4) and derivate: 
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The second term on the right member can be further developed as follows: 

where we have integrated by parts and used Eq. (3). Then, it is possible to rewrite (26): 
On the other hand, the derivation of the general expression for the Hartree energy leads to: 
Equating (28) and (29) we find the following relation: 

Ultimately, replacing (30) into (28) we obtain the final result, 

dEH 1 [ deini) ^ [ A.( \^Piotij^) ^ ^Q1^ 
mrr-^l^^^^'^^'^^ ^^ + y'^W^^^r. (31) 
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TABLE I: Solvation free energies (kcal/mol) for selected molecules and ions in water, calculated 
with this model using a dielectric function e = e{p) determined by the electron density,'^ and 
e = e(7(Ri)) determined by the atomic positions with /?=2.4 (see text). Experimental values'' and 
results from PCM^ are also shown. 



Expt. e = e(/5) e = e(7) PCM 

H2O -6.3 -8.4 -8.8 -5.4 

CH3CONH2 -9.7 -10.5 -8.0 -4.6 

CH3NH+ -73 -81.0 -81.9 -65.1 

NO3 -65 -57.8 -60.6 -62.6 

Cr -75 -66.9 -68.6°^ -72.6 



From reference^. ^ From referencea^i^. Obtained with the Polarizable Continuum Model as 
implemented in Gaussian 03.— i2I The van der Waals radius for ionic chlorine was set equal to 

2.059 A, from reference^. 
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TABLE II: Adsorption energies for water on the anatase (101) surface in the gas phase and in 
solution, in both the molecular and dissociative configurations at different coverages. Values are 
given in kcal per H2O molecule. 



Molecular Dissociative 

e = 0.25 e = l 9 = 0.25 9 = 1 

Gas phase -19.4 -17.8 -1.3 -7.4 

Solution - -3.0 - 6.6 
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FIG. 1: Total energy during two different Car-Parrinello electronic minimizations in a continuum 
solvent (the structure is the Ti02 slab described in section III). The squares indicate the results for 
a dielectric function depending on the charge density, whereas the circles correspond to the same 
computational experiment, but removing from the total potential (see text). Both runs used 
the same calculation parameters and were restarted from the wavefunction converged in vacuum. 
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FIG. 2: Variation of and of the effective potential along the z-direction (perpendicular to the 
slab) at a selected point on the x-y plane, corresponding to the position of one of the exposed Ti 
atoms of the surface. The dotted line represents the permittivity, referenced on the right axis. The 
arrows point to the approximate location of each of the titanium layers. 
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FIG. 3: The permittivity around an oxygen atom as a function of the distance, for different values 
of the parameter /?, according to the position-dependent dielectric defined in Eqs. (9) and (10). 
The transition, centered around the van der Waals radius, becomes sharper as (3 is increased. 




FIG. 4: Contour plot of the dielectric function e(r) in a supercell containing a four layers slab 
representing the anatase (101) face of Ti02. An isosurface corresponding to e(r)=1.4 is displayed 
in yellow. In (a) the solvent percolates through the surface, whereas in (b) it is excluded from the 
slab by virtue of the artifact given in relation (19). 
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FIG. 5: Interatomic distances Oa-Ha and 0;,-Ha during a molecular dynamics simulation of a proton 
transfer process, from an adsorbed water molecule to an hydroxide ion in a continuum solvent. 
Oa and O5 are the oxygen atoms of the initially absorbed water molecule and the hydroxide, 
respectively, and Hq is the exchanged proton. At the beginning of the simulation the OH^ group 
is H-bonded to the water molecule. The geometry of the system is shown in the inset for selected 



time steps. 
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FIG. 6: Interatomic distance 0{,-Ha during the molecular dynamics simulations of an adsorbed 
water molecule in the presence of an hydroxide ion initially situated 4A above the surface. 
is the oxygen atom of the hydroxide and is a proton of the water molecule, which is rapidly 
exchanged during the gas phase dynamics. The starting geometry is shown on the left, while the 
upper and lower figures on the right depict the atomic structures after 1 ps of dynamics in vacuum 
and in solution, respectively. 
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FIG. 7: Electron density integrated on the xy plane, and displayed as a function of the z coordinate, 
for the hydroxyl anion situated 4A above the surface. The approximate positions of the upper Ti 
layer and of the 0H~ ion are indicated with arrows. Enlargement of the interfacial region shows a 
depletion of the electron density in solution with respect to vacuum. There is an average difference 
of about 0.5 e between the two curves, owing to the stabilization of the electronic charge of the 
hydroxyl anion surrounded by the dielectric. 
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